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Universality is a well-established central concept of equilibrium physics. However, in systems 
far away from equilibrium a deeper understanding of its underlying principles is still lacking. Up 
to now, a few classes have been identified. Besides the diffusive universality class with dynamical 
exponent z = 2 another prominent example is the superdiffusive Kardar-Parisi-Zhang (KPZ) class 
with 2 = 3/2. It appears e.g. in low-dimensional dynamical phenomena far from thermal equilibrium 
which exhibit some conservation law. Here we show that both classes are only part of an infinite 
discrete family of non-equilibrium universality classes. Remarkably their dynamical exponents Za 
are given by ratios of neighbouring Fibonacci numbers, starting with either zi = 3/2 (if a KPZ 
mode exist) or 2 i = 2 (if a diffusive mode is present). If neither a diffusive nor a KPZ mode 
are present, all dynamical modes have the Golden Mean z = (1 + ■\/5)/2 as dynamical exponent. 
The universal scaling functions of these Fibonacci modes are asymmetric Levy distributions which 
are completely fixed by the macroscopic current-density relation and compressibility matrix of the 
system and hence accessible to experimental measurement. 


The Golden Mean (/? = 1/2 -|- •\/5/2 Ri 1.61803.., also 
called Divine Proportion, has been an inspiring number 
for many centuries. It is widespread in nature, i.e. ar¬ 
rangements of petals of the flowers and seeds in the sun¬ 
flower follow the golden rule [I|. Being considered an 
ideal proportion, the Golden Mean appears in famous ar¬ 
chitectural ensembles such as Parthenon in Greece, Giza 
Great Pyramids in Egypt, or Notre Dame de Paris in 
France. Ideal proportions of the human body follow the 
Golden Rule. 

Mathematically, the beauty of the Golden Mean num¬ 
ber is expressed in its continued fraction representation: 
All the coefficients in the representation are equal to 
unity, 


Systematic truncation of the above contin¬ 
ued fraction gives the so-called Kepler ratios, 
1/1, 2/1, 3/2, 5/3,8/5,..., which approximate the Golden 
Mean. Subsets of denominators (or numerators) of the 
Kepler ratios form the celebrated Fibonacci numbers, 
Fi = 1,1, 2, 3, 5, 8 ,.., such that Kepler ratios are ratios 
of two neighbouring Fibonacci numbers. As well as the 
Golden Mean, Fibonacci ratios and Fibonacci numbers 
are widespread in nature [l|. 

The occurrence of the Golden Mean is not only in¬ 
teresting for aesthetic reasons, but often indicates the 
existence of some fundamental underlying structure or 


symmetry. Here we demonstrate that the Divine Pro¬ 
portion, as well as all the truncations (Kepler ratios) of 
the continued fraction appear as universal numbers, 
viz., the dynamical exponents, in low-dimensional dy¬ 
namical phenomena far from thermal equilibrium. The 
two well-known paradigmatic universality classes, Gaus¬ 
sian diffusion with dynamical exponent z = 2 d, d and 
the Kardar-Parisi-Zhang (KPZ) universality class with 
0 = 3/2 enter the Kepler ratios hierarchy as the first 
two members of the family. 

The universal dynamical exponents in the present con¬ 
text characterize the self-similar space-time fluctuations 
of locally conserved quantities, characterizing e.g. mass, 
momentum or thermal transport in one-dimensional sys¬ 
tems far from thermal equilibrium Q ■ The theory of non¬ 
linear fluctuating hydrodynamics (NLFH), has recently 
emerged as a powerful and versatile tool to study space- 
time fluctuations, and specifically the dynamical struc¬ 
ture function which describes the behaviour of the slow 
relaxation modes, and from which the dynamical expo¬ 
nents can be extracted Q. 

The KPZ universality class m has been shown to 
explain the dynamical exponent observed in interface 
growth processes as diverse as the propagation of flame 
fronts [3, d , the growth of bacterial colonies dl ; or the 
time evolution of droplet shapes such as coffee stains 0 
where the Gaussian theory fails. The dynamical struc¬ 
ture function originating from the one-dimensional KPZ 
equation has a non-trivial scaling function obtained ex¬ 
actly by Prahofer and Spohn from the totally asymmetric 
simple exclusion pro cess (TASEP) and the polynuclear 
growth model [l5l.H6l| and was beautifully observed in ex- 
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periments on turbulent liquid crystals [i3[i3- Theoret¬ 
ical treatment, both numerical and analyticahof generic 
model systems with Hamiltonian dynamics [l9|, anhar- 
monic chains 1^. 1^ and lattice models for driven dif¬ 
fusive systems , have demonstrated an extraordi¬ 

nary robust universality of fluctuations of the conserved 
slow modes in one-dimensional systems. 

Despite this apparent ubiquity, dynamical exponents 
different from z = 2orz = 3/2 were observed frequently. 
Usually it is not clear whether this corresponds to gen¬ 
uinely different dynamical critical behaviour or is just a 
consequence of imperfections in the experimental setting. 
Moreover, recently a new universality class with dynam¬ 
ical exponent z = 5/3 for the heat mode in Hamiltonian 
dynamics [l^ was discovered, followed by the discovery 
of some more universality classes in anharmonic chains 
2C, 1^ and lattice models for driven diffusive systems 
2^ [^ . What is lacking, even in the conceptually sim¬ 
plest case of the effectively one-dimensional systems that 
we are considering, is the understanding of the plethora of 
dynamical non-equilibrium universality classes within a 
larger framework. Such a framework exists e.g. for two- 
dimensional critical phenomena in equilibrium systems 
where the spatial symmetry of conformal invariance to¬ 
gether with internal symmetries give rise to discrete fam¬ 
ilies of universality classes in which all critical exponents 
are simple rational numbers. 

It is the aim of this article to demonstrate that dis¬ 
crete families of universality classes with fractional crit¬ 
ical exponents appear also far from thermal equilib¬ 
rium. This turns out to be a hidden feature of the 
NLFH equations which we extract using mode-coupling 
theory. It is remarkable that one finds dynamical ex¬ 
ponents Zq, which are ratios of neighbouring Fibonacci 
numbers {1,1, 2, 3, 5,8,...} defined recursively as = 
Fn-i -f Fn- 2 - The first two members of this family are 
diffusion (z = 2 = F 3 /F 2 ) and KPZ (z = 3/2 = F 4 /F 3 ). 
The corresponding universal scaling functions are com¬ 
puted and shown to be (in general asymmetric) Za-stable 
Levy distributions with parameters that can be com¬ 
puted from the macroscopic current-density relation and 
compressibility matrix of the corresponding physical sys¬ 
tem and which thus can be obtained from experiments 
without detailed knowledge of the microscopic proper¬ 
ties of the system. The theoretical predictions, obtained 
by mode coupling theory, are confirmed by Monte-Carlo 
simulations of a three-lane asymmetric simple exclusion 
process which is a model of driven diffusive transport of 
three conserved particle species. 


I. NONLINEAR FLUCTUATING 
HYDRODYNAMICS 

We consider a rather general interacting non¬ 
equilibrium system of length L described macroscopically 
by n conserved order parameters p\(x, t) with stationary 
values p\ and associated macroscopic stationary currents 


j\{pi,..., Pn) and compressibility matrix K with ma¬ 
trix elements Kxp. = ^ ((^a - PxL){Nfj, - p^L)) where 

N\ = dxpx(x,t) are the time-independent conserved 
quantities. 

The starting point for investigating density fluctua¬ 
tions u\{x,t) := px{x,t) — px in the non-equilibrium 
steady state are the NLFH equations Q 

dtu =-dj; ^Ju + ^{u\H\u) - d^Du + (2) 

where J is the current Jacobian with matrix elements 
= djx/dpfj,, iL is a column vector whose en¬ 
tries = F[^ are the Hessians with matrix el¬ 
ements = d'^jx/{dp^dpv) and the bra-ket nota¬ 

tion represents the inner product in component space 
{u\{H)x\u) = with (u | = 

\u) = u. The diffusion matrix I? is a phenomenological 
quantity. The noise term does not appear explicitly 
below, but plays an indirect role in the mode-coupling 
analysis. The product JK of the Jacobian with the com¬ 
pressibility matrix K is symmetric whic h g uarantees 
a hyperbolic system of conservation laws [^. We ig¬ 
nore possible log arithmic corrections arising from cubic 
contributions |26l |. 

This system of coupled noisy Burgers equations is con¬ 
veniently treated in terms of normal modes (j) = Ru 
where RJR~^ = diag(uQ,) and the transformation ma¬ 
trix R is normalized such that RKR^ = I. The 
eigenvalues Va of J are the characteristic velocities of 
the system. From m one thus arrives at dt4>ct = 

-da: (va(l)a + - da:iD(f)a + withD = 

RDR~^, B = RB and the mode coupling matrices 

G“ = i ^ R^piR-^fH^R-^ (3) 

whose matrix elements we denote by G'p^. 


II. COMPUTATION OF THE DYNAMICAL 
STRUCTURE FUNCTION 

The dynamical structure function describes the sta¬ 
tionary fluctuations of the conserved slow modes and 
is thus a key ingredient for understanding the interplay 
of noise and non-linearity and their role for transport 
far from equilibrium. We focus on the case of strict 
hyperbolicity where all Va are pairwise different and 
study the large scale behaviour of the dynamical struc¬ 
ture function S°‘^{x,t) = {(j)a{x,t)(j)p{0,0)). Since all 
modes have different velocities only the diagonal elements 
Sa{x,t) := S°‘°‘{x,t) are non-zero for large times. Mode 
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coupling theory yields @ 

dtSa{x,t) = {-Vad^ + Dadl)Sa{x,t) 

+ [ ds [ dySa{x - y,t - s)dlMaaiy,s) (4) 
Jo Jr 

with the diagonal element Da '■= Daa of the phenomeno¬ 
logical diffusion matrix for the eigenmodes and the mem¬ 
ory kernel Maa{y,s) = ‘ 2 -Y.i 3 ,-yi.G‘^~ifSfi{y,s)S^[y,s). 
The task therefore is to extract for arbitrary n the large¬ 
time and large-distance behaviour from this non-linear 
integro-differential equation. 

Remarkably, these equations can be solved exactly in 
the long-wavelength limit and for t —> oo by Fourier and 
Laplace transformation (see Appedices). Using a suit¬ 
able scaling ansatz for the transformed structure func¬ 
tion then allows to analyze the small-p behaviour from 
which the dynamical exponents can be determined. We 
find that different conditions arise depending on which 
diagonal elements of the mode-coupling matrices vanish. 


III. THE FIBONACCI FAMILY OF 
DYNAMICAL UNIVERSALITY CLASSES 

A. Fibonacci case 


First, we consider the case where the self-coupling 
is nonzero for one mode only, e.g. G\i ^ 0. For all 
other modes a > 1 we assume a single nonzero coupling 
to the previous mode, so Ga-i,a-i Oj — 

0 for fj ^ a — 1. Then, as follows from our analysis 
(see Appendices), we find the following recursion for the 
dynamical exponents: 

1 

Za — 1 H- (5) 

^OL — 1 


with zi = 3/2. 

The dynamical structure function in momentum space 
is proportional to the Zo,-stable Levy distribution with 
maximal asymmetry (t“ = ±1, see and eq. (IA4I) be¬ 
low. The sign of the asymmetry depends whether the 
mode (a—1) has bigger or smaller velocity than the mode 
a, (t“ = —sgn(uQ, — Va-i)- The dynamical exponents ([5]) 
form a sequence of rational numbers 


^a. 


Fg+O 

Fa+2 


( 6 ) 


which are consecutive ratios of neighbouring Fibonacci 
numbers Fq, defined by Fa = Fa-i + Fa -2 with initial 
values Fq = 0, Fi = 1, which converge exponentially 
to the Golden Mean (p := i(l -I- v^) r; 1.618, as first 
observed by Kepler in 1611 in a treatise on snow flakes. In 
a model with n conservation laws, one has the Fibonacci 
modes with dynamical exponents {3/2, 5/3, 8/5,..., z„}. 

Finally, we remark that if mode 1 is diffusive rather 
than KPZ, then we find the same sequence ([5]) of expo¬ 
nents except that it starts with zi = F 2 = 2. 



FIG. 1: The scaling functions (bottom) and dynamical expo¬ 
nents are related to the structure of the mode coupling ma¬ 
trices G“ (top). The table shows the dynamical exponents Za 
in the case n = 2, see eq. dMI- The symbols * and * denote 
non-zero elements. Red symbols correspond to self-coupling, 
black symbols to couplings to other modes. Matrix elements 
not indicated can take any value. The colors in the table cor¬ 
respond to the colors of the graphs of the scaling functions. 


In Fig. [T]we show some representative examples of the 
scaling functions which are quite different in shape. Fur¬ 
thermore the relation between the exponents Za , deter¬ 
mined by eq. (IA2I) . and the mode coupling matrices G“ 
is illustrated for the case n = 2. 


B. Golden Mean case 

As second representative example we consider the case 
where all self-coupling coefficients vanish, G“q, = 0 for 
all a, while each mode has at least one nonzero cou¬ 
pling to another mode, G^^ yf 0 for some jJ ^ a. Then, 
(IS|) reduces to Za = \ F^jzp for all modes a, fJ. The 
unique solution of this equation is the Golden Mean 
-Za = = (1 + ■\/5)/2 for all a. The scaling functions 
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FIG. 2: Schematic drawing of three particle species drifting 
inside a nanotube. Due to the interaction between the parti¬ 
cles and with the walls one expects a non-linear current den¬ 
sity relation. 


(see Supporting Information section) are proportional to 
(p-stable Levy distributions with parameters fixed by the 
collective velocities and the mode-coupling coefficients. 
The asymmetry of the fastest right-moving (left-moving) 
mode is predicted to be /3 = — 1 (/3 = 1). 


IV. SIMULATION RESULTS 

To check the theoretical predictions for the two cases 
we simulate mass transport with three conservation laws, 
i.e., three distinct species of particles. To maintain a far- 
from-equilibrium situation a driving force is applied that 
leads to a constant drift superimposed on undirected dif¬ 
fusive motion. This is a natural setting for transport of 
charged particles in nanotubes, see Fig. [2] for an illustra¬ 
tion, where a direct measurement of the stationary parti¬ 
cle currents is experimentally possible [2^. However, due 
to the universal applicability of NLFH the actual details 
of the interaction of the particles with their environment 
and the driving field are irrelevant for the theoretical de¬ 
scription of the large-scale dynamics. Hence for good 
statistics we simulate a lattice model for transport which 
represents a minimal realization of the essential ingredi¬ 
ents, namely a non-linear current-density relation for all 
three conserved masses. 

Our model is the three-species version of the multi-lane 
totally asymmetric simple exclusion process [2^. Parti¬ 
cles hop randomly in field direction on three lanes to their 
neighbouring sites on a periodic lattice of 3 x L sites with 
rates that depend on the nearest-neighbour sites. Lane 
changes are not allowed so that the total number of par¬ 
ticles on each lane is conserved. Due to excluded-volume 
interaction each lattice site can be occupied by at most 
one particle. Thus the occupation numbers of site k 
on lane A take only values 0 or 1. The hopping rate 



FIG. 3: Space-time propagation of three normal modes in 
the three-lane model. The modes (from left to right) are the 
Fibonacci mode with 2 = 8/5 (mode 3), the KPZ mode with 
z = 3/2 (mode 1), and the Fibonacci mode with z = 5/3 
(mode 2). The physical and simulation parameters are given 
in the Appendices. 


from site k on lane A to site fc -I- 1 on the same lane is 
given by 

^ + ”i+i) 

fri: m/A} 

with a species dependent drift parameter b\ and symmet¬ 
ric interaction constants = 7 ^a- Hopping attempts 
onto occupied sites are rejected. The conserved quanti¬ 
ties are the three total numbers of particles Nx on each 
lane with corresponding densities px = Nx/L. 

The stationary distribution of our model factorizes 
[ 2 ^ and thus allows for the exact computation of the 
macroscopic current-density relations jxipi^ P 2 , Ps) and 
the compressibility matrix K{pi, p 2 , pz). Furthermore, 
because there is no particle exchange between lanes, the 
compressibility matrix is diagonal with elements denoted 
by Kx- One has 

JA = Pa(I-Pa) Ua+ X! (8) 

V / 

Kx = px{l-px)- (9) 

The diagonalization matrix R and the mode-coupling 
matrices G“ are fully determined by these quantities. 

According to mode-coupling theory three different 
Fibonacci-modes with zi = 3/2, Z 2 = 5/3, 2:3 = 8/5 
occur e.g. when G\i ^ 0 , ^ 0 , G 22 ^ 0 , and 

G 22 = GI 3 = G 33 = 0. For our simulation we com¬ 
pute numerically densities, bare hopping rates and in¬ 
teraction parameters to satisfy these properties as de¬ 
scribed in the Appendices. For this choice the velocities 
of the normal modes are vi = 0.592315, V 2 = 0.0281578, 
U 3 = 1.58226 which ensures a good spatial separation 
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FIG. 4: Left Panel: Vertical least squares fit of the numeri¬ 
cally obtained dynamical structure function for the Fibonacci 
8/5-mode (points), at time t = 1000 with a 8/5-stable Levy 
distribution, maximal asymmetry —1 and theoretical center 
of mass (line), predicted by the mode coupling theory. The 
only fit parameter is the scale parameter of the Levy sta¬ 
ble distribution. The simulation results agree very well with 
the asymptotic theoretical result already for moderate times. 
Right Panel: Insets show closeups of the peak region and tail 
regions, according to a colour code. Every 10-th datapoint is 
plotted to improve the visibility of the data. The statistical 
error 699 % with 99% confidence bound is for every data point 
smaller than 1.6299 • 10“®. 

after quite small times. The propagation of the three 
normal modes (Fig. [3]) with the predicted velocities is 
observed with an error of less than 10“^. Moreover, the 
numerically obtained dynamical structure function for 
mode 3 shows a startling agreement with the theoreti¬ 
cally predicted Levy scaling function with z = 8/5 and 
maximal asymmetry, see Fig. 21 It takes longer for the 
other two modes (KPZ mode and Levy stable 5/3 mode) 
to reach their asymptotic form, which we argue is due to 
the much smaller respective couplings, (Gii/G 22 )^ ^ 
(G?i/Gi2)2 « 1. 

In order to observe three Golden Mean modes it is 
sufhcient to require that each mode has zero self-coupling 
and at least one nonzero coupling to other modes. This 
can be achieved with the set of parameters given in 
Appendices which lead to the velocities vi = 1.83149, 
V 2 = 0.762688, V 3 = 0.326778 of the normal modes. 
The propagation of the three normal modes with the 
predicted velocities is observed, approaching for large 
times a very small relative error of about 10“^. The 
structure function for the fastest mode 1 converges to its 
asymptotic shape faster than for the other modes, due 
to the large coupling coefficient G 33 . In Fig. [S] we show 
a scaling plot of the measured structure function for 
mode 1 with dynamical exponent z = tp = (l + v^) /2 
together with a fitted to a (^-stable Levy function 
(lAl) with maximal asymmetry /3 = —1 as predicted 
by the theory. The data collapse shows a striking 
agreement between the measured and theoretical scaling 


FIG. 5: Scaling plot of the measured structure function of 
mode 1 with dynamical exponent z = (p = (l -I- %/5) /2 for the 
Golden Mean case, htted to a :^-stable Levy distribution with 
maximal asymmetry —1 (see Eq. dH). The scale parameter 
El for the Levy-stable distribution and the center of mass ve¬ 
locity vi are obtained by a vertical least square fit. Fitted pa¬ 
rameters are vifit = 1.83107T0.00009 and = l.lTO.Ol. 
The fitted velocity differs by 0.02% from the theoretical 
velocity. 

function. Alternatively, the dynamical exponent Za 
can be derived from the maximum of the structure 
function, which scales as max(S'i(a:, t)) = const • . 

We obtain z si 1.63 which differs from the predicted 
value z = tp = {1 + •\/5)/2 by less than 0.8%. 


V. DISCUSSION 

Our work demonstrates that non-equilibrium phenom¬ 
ena are much richer than just the diffusive and KPZ 
universality suggest. We have established that in non¬ 
equilibrium phenomena governed by non-linear fluctuat¬ 
ing hydrodynamics with n conservation laws mode cou¬ 
pling theory predicts a family of dynamical universality 
classes with dynamical exponents given by the sequence 
of consecutive Kepler ratios ([ 6 |) of Fibonacci numbers. 
With slightly modified initial conditions on G}^ this re¬ 
sult is easily generalized for the case when the first mode 
a = 1 is diffusive. Then the sequence of dynamical ex¬ 
ponents becomes shifted by one unit with respect to ([HI) . 
On the other hand, if all self-couplings vanish, but at 
least one other diagonal element G^^ of the mode cou¬ 
pling matrix is non-zero,one has as unique solution for 
all modes a the fixed point value Za = Zoo = which is 
the Golden Mean. 

For general mode coupling matrices all critical expo¬ 
nents can be computed (from (IA2I) in Appendices). The 
scaling functions of the non-diffusive and non-KPZ modes 
are asymmetric Levy distributions whose parameters 
are completely determined by the macroscopic current- 
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density relation and compressibility matrix of the system. 

For 1+1 dimensional systems out of equilibrium this is 
the first time, to our knowledge, that an infinite family of 
discrete universality classes is found. Recalling that 1+1 
dimensional non-equilibrium systems with short-range 
interactions can be mapped onto two-dimensional equi¬ 
librium systems (with the time evolution operator play¬ 
ing the role of the transfer matrix) one is reminded of the 
discrete families of conformally invariant critical equilib¬ 
rium systems in two space dimensions [13, IMj. We do 
not know whether there is any mathematical link, but 
the analogy is suggestive in so far as conformal invari¬ 
ance is a local symmetry of spatially isotropic systems 
with z = 1 (which happens to be the lowest order Kepler 
ratio) while z > 1 corresponds to strongly anisotropic 
systems for which also local symmetry groups are known 
to exist [33| . 

Since an infinite number of lanes of coupled one¬ 
dimensional systems corresponds to a two-dimensional 
system, it is intriguing to observe that the Golden Mean 
is close to the numerical value z = 1.612 —1.618 of the dy¬ 
namical exponent of the 2 + 1-dimensional KPZ-equation 
[ 33,113 • The scaling function of the 2 + 1-dimensional 
KPZ-equation, however, is not Levy [33| . 

In order to observe and distinguish between the differ¬ 
ent new classes highly precise experimental data will be 
required. E.g. in the Fibonacci case the dynamical expo¬ 
nents converge quickly to the Golden Mean. A feature 
which might be easier to observe experimentally is the 
scaling function itself, which for higher Fibonacci ratios 
5/3, 8/5,... usually has a strong asymmetry (see Figs. [H 
|4]and[5]) while KPZ and Gauss scaling functions are sym¬ 
metric. Growth processes which can be mapped on exclu¬ 
sion processes with several conservation laws, might be 
potentially suitable candidates for an experimental veri¬ 
fication, see e.g. HO for an example of a system with 
one conservation law. 


Appendix A: Computation of the dynamical 
structure function 

The mode coupling equations O can be solved in 
the scaling limit by applying a Fourier transform (FT) 
f{x) —>■ f{p) and a Laplace transform (LT) /(f) —>• /(w). 
For more details we refer to [131 where the case n = 2 
of two conservation laws has been treated. After making 
the scaling ansatz 

Saip^Cboi) —p 9a{Ca') (-^1) 

for the transformed dynamical structure function where 
Sa{p, 0) = l/v^ and (a = WqIpI we are in a position 
to analyze the small-p behaviour. One has to search for 
dynamical exponents for which the limit p —>■ 0 is non¬ 
trivial, which requires a self-consistent treatment of all 
modes. We find that different conditions arise depending 
on which diagonal elements of the mode-coupling ma¬ 
trices vanish. To characterize the possible scenarios we 


define the set !„ := {/3 : ^ 0} of non-zero diagonal 

mode coupling coefficients. Through power counting one 
obtains 


Zqc 


2 

3/2 

min/3gi^ 



if la = 0 

if a S la 

else 


and the domain 


I < Za<2 Va. 


(A2) 


(A3) 


for the possible dynamical exponents. 


In the Fibonacci case, the dynamical structure function 
of mode a in momentum space has the scaling form 

S^{p^t) = (A4) 

V 27r 

with inverse time scales Ea- The dynamical exponents 
then satisfy the recursion ([5]). Up to the normalization 
Xj^ph^ the scaling form (IA4I) is a a-stable Levy distri¬ 
bution (27j| . 


Appendix B: Simulation algorithm 

For the Monte-Garlo simulation of the model we choose 
a large system size L > 5 • 10® which avoids finite-size ef¬ 
fects. At time t = 0, N\ particles are placed on each lane 
according to the desired initial state. One Monte-Garlo 
time unit consists of 3 • L • r* random sequential update 
steps where r* = max{rfc('^^}: In each update step a 
bond + 1) is chosen randomly with uniform 

distribution. If = 1 then the particle at 

site k is moved to A: + 1 with probability /r* where 

r* is the maximal value that the can take among 
all possible particle configurations on the neighbouring 
lanes. If ^^"^^(1 — = 0 the particle configuration 

remains unchanged. 


Appendix C: Simulation of the dynamical structure 
function 

In order to determine the dynamical structure function 
we initialize the system by placing Nx particles uniformly 
on each lane A. This yields a random initial distribution 
drawn from the stationary distribution of the process. 
No relaxation is required. 

Then we use translation invariance and compute the 
space- and time average 

M L 

j=i 

(Cl) 
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To avoid noisy data of we take in m the system 
size L and the time average parameter M sufficiently 
large. In order to obtain we average over P in¬ 

dependently generated and propagated initial configura¬ 
tions of The error estimates for are cal¬ 

culated from the P independent measurements. From 
we compute the structure function of the normal 
modes by transformation with the diagonalizing matrix 
R determined by dl]) and dl]). 

To obtain model parameters for three different 
Fibonacci-modes with zi = 3 / 2,22 = 5 / 3,23 = 8/5 we 
solve the equations given in the text after Eq. numer¬ 
ically with a C-program performing direct minimization 
of the absolute values of the targeted G-elements until 
the given tolerance value (10“®) is reached. The data 
shown here for the three mode case have been obtained 
from simulations with densities pi = 0.2, p 2 = 0.25, pi = 
0.3, bare hopping rates bi = 0.613185, 62 = 0.425714, 
63 = 0.799831 and interaction parameters 712 = 1.36145, 
723 = 3.69786, 713 = 0.143082 for which the needed re¬ 
lations are satisfied. This choice of parameters yields 
Gji = 0.322507, = -0.15, = 1.04547, while 

the the absolute values of G 22 , G§ 3 , G 33 are smaller than 
10“®. Besides these physical parameters, the simula¬ 
tion parameters for the Fibonacci modes (Fig. [3l |4]) are 
L = 5 • 10^ T = 250, M = 1400, P = 98. 

For the Golden Mean case (Fig. [S]) the set of param¬ 
eters Pi = 0.2, p 2 = p 3 = 0.25, 712 = 0.0082334758646, 
723 = 1.68447706968, 713 = 3.72140740146, and h = 
0.905073261248, 62 = 0.86, and 63 = 1.18875738638. 
This leads to G ^2 = 0.405702, G^g = 0.929315, Gf^ = 
-0.104141, Gi 3 = -0.208477, Gfi = -0.182467, GI 2 = 
0.271246, while the absolute values of Gli,G 22 ,Gl^ are 
smaller than 10“®. The simulation parameters for the 
Golden Mean case are L = 5 • 10®, t = 750, M = 30 and 
P = 303. 


Supporting Information 

Remarkably, the mode coupling equations eq. o can 
be solved exactly in the scaling limit by Fourier and 
Laplace transformation. To this end we define the Fourier 
transform (FT) as 


f{p) := ^ r dxe-^pyix), (G2) 

V J —00 

and the Laplace transform (LT) as 

poo 

/» := / dte--‘/(0- (C3) 

Jo 

With Da (p) = iVaP + DaP^ we obtain from eq. (4) of the 
paper in momentum-frequency space 


5'a(p,w) 


_ Sa{p,d) _ 

W -I- Da{p) + CaaiP,(^) 


(G4) 


with memory kernel 


G„„(p,a;) = 2^(G|^)V 

/ 3,7 

and Sa{p,0) = l/\/27r. 

Next we introduce Wq, := 
ansatz 



dse"'^'* 


dqSp{q,s)Sj{p-q, s). 
(G5) 


ujpiVaP and make the scaling 


Sa{p,0Ja)=P "'“5a(Ca) (C6) 

with Co, = a)o|p|~^“- Having in mind systems with short- 
range interactions we anticipate that all modes spread 
subballistically, i.e., Za > I for all a. Using strict hyper- 
bolicity one obtains after some substitutions of variables 


gaiCa) = lim 

p —>-0 




/3/a 


(C7) 


with := \va — U/3|sgn[p(uQ — vp)] and 

Qap = 2{G^pfT (^1 - ^j > 0. (G8) 

where 

/ OO 

dpHp)f{-p)- (C9) 

-00 

With cr“^ = sgn[p(uo — vp)] one has 

^ = sin |ua - 




(CIO) 


Now we are in a position to analyze the small-p be¬ 
haviour. One has to search for dynamical exponents 
for which the limit p —>■ 0 is non-trivial, which is de¬ 
termined by the smallest power of p in (IC7I) . This has 
to be done self-consistently for all modes. We find that 
different conditions arise depending on which diagonal 
elements of the mode-coupling matrices vanish. In order 
to characterize the possible scenarios we define the set 
Iq, := {/3 : G'pp ^ 0 } of non-zero diagonal mode cou¬ 
pling coefficients. One obtains from (ICTl) through power 
counting the system of equations 




2 

3/2 

min;3Gic 



if la = 0 

if a £ la 

else 


and the domain 


(Cll) 


1 < Za < 2 Va. (C12) 

for the possible dynamical exponents. 
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